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g ■ ABSTRACT 

O 



By numerically integrating the orbits of the giant planets and of test particles 
^ ! over a period of four billion years, we follow the evolution of the location of the 

midplane of the Kuiper belt. The Classical Kuiper belt conforms to a warped 
sheet that precesses with a 1.9 Myr period. The present-day location of the 
^ ! Kuiper belt plane can be computed using linear secular perturbation theory: the 

local normal to the plane is given by the theory's forced inclination vector, which 
is specific to every semi-major axis. The Kuiper belt plane does not coincide with 
^ [ the invariable plane, but deviates from it by up to a few degrees in stable zones. 

For example, at a semimajor axis of 38 AU, the local Kuiper belt plane has an 
^ [ inclination of 1.9 degrees and a longitude of ascending node of 149.9 degrees when 

K* I referred to the mean ecliptic and equinox of J2000. At a semimajor axis of 43 

QQ I AU, the local plane has an inclination of 1.9 degrees and a nodal longitude of 78.3 

^ I degrees. Only at infinite semimajor axis does the Kuiper belt plane merge with 

^ I the invariable plane, whose inclination is 1.6 degrees and nodal longitude is 107.7 

O ' degrees. A Classical Kuiper belt object keeps its inclination relative to the Kuiper 

00 ' 

Q I belt plane nearly constant, even while the latter plane departs from the trajectory 

predicted by linear theory. The constancy of relative inclination reflects the 

undamped amplitude of free oscillation; that is, the homogeneous solution to the 



H ' forced harmonic oscillator equation retains constant amplitude, even while the 

inhomogeneous solution cannot be written down accurately because the planetary 
forcing terms are chaotic. Current observations of Classical Kuiper belt objects 
are consistent with the plane being warped by the giant planets alone, but the 
sample size will need to increase by a few times before confirmation exceeds 3a 
in confidence. In principle, differences between the theoretically expected plane 
and the observed plane could be used to infer as yet unseen masses orbiting the 
Sun, but carrying out such a program would be challenging. 
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Subject headings: comets: general — Kuiper belt — solar system: general 
celestial mechanics 



1. INTRODUCTION 

If we could map, at fixed time, the instantaneous locations in three-dimensional space of 
all Kuiper belt objects (KBOs), on what two-dimensional surface would the density of KBOs 
be greatest? We call this surface the plane of the Kuiper belt (KBP), though by "plane" we 
do not mean to imply that the KBP is flat (we shall find that it is not). The KBP depends 
on the mass distribution of the solar system — principally, the orbits and masses of the giant 
planets. There are as many different KBPs as there are dynamical classes of KBO, since 
each class of object feels a distinct time-averaged force. Here we study the KBP defined by 
Classical KBOs: objects whose fairly circular, low inclination orbits are not in any strong 
mean-motion resonance with Neptune (see Elliot et al. 2005 for a classification scheme). 

In principle, theoretical determination of the KBP would help observers to discover new 
KBOs. Conversely, by measuring differences between the theoretical KBP and the actual 
KBP, we might hope to infer the presence of solar system bodies as yet undetected ( "Planet 
X"; see Gaudi & Bloom 2005 for a summary of current limits). 

There is disagreement regarding the location of the KBP. Brown & Pan (2004) analyzed 
the instantaneous proper motion vectors of hundreds of KBOs irrespective of dynamical class 
and concluded, with greater than 3a confidence, that the KBP did not coincide with the 
invariable plane (IP, the plane perpendicular to the total angular momentum vector of the 
solar system). They argued that the observed KBP was consistent instead with the forced 
plane given by linear secular perturbation theory. We will refer to this plane as the BvWP, 
after Brouwer & van Woerkom (1950), who developed a linear secular theory for the motions 
of all eight of the major planets. Their theory, in turn, has its origin in the Laplace-Lagrange 
equations (see, e.g., Murray Sz Dermott 1999). The BvWP is a warped and time-variable 
surface whose properties we review in §5J 

By contrast, Elliot et al. (2005) found that the plane determined by Classical KBOs 
that were observed over multiple epochs was more consistent with the IP (< la difference) 
than with the BvWP (< 2-3cr difference). They listed some arguments, none conclusive, for 
why the IP might be preferred over the BvWP. The low order of the BvW theory, and its 
inability to account for time variations in semi-major axes, are causes for concern. 

We seek to resolve this disagreement using numerical orbit integrations. In §2l we review 
the linear theory and how it equates the KBP with the BvWP. In §31 numerical integrations 
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lasting the age of the solar system are used to reveal the theoretical location of the KBP. 
In mi we compare theory against current observations of the KBP. A summary is given in 
§3 including a brief comment on the prospects for detecting an unseen, outer solar system 
planet using the KBP. 



2. LINEAR SECULAR THEORY 

We study the secular evolution of Classical KBOs by solving the Laplace-Lagrange 
equations of motion, which neglect all terms higher than second order in orbital eccentricity 
(e) and orbital inclination [i). The solution is detailed in the textbook by Murray & Dermott 
(1999); we provide a summary here. 

We start by describing the motions of the planets. Define an inclination vector i 
= {q,p) = {icosQ,ismQ), where i and Q equal the inclination and longitude of ascending 
node. Lagrange's equations governing the inclination vector for the jth planet are 

1 dR, 



^^ - . .2 



Ujaj dpj 

P=^^, (2) 

■^ rijOj dqj 

where rij, aj, and Rj are, respectively, the mean motion, semi-major axis, and disturbing 
function. We consider only the four giant planets so that j = 1,2,3,4 represents Jupiter, 
Saturn, Uranus, and Neptune, respectively. The disturbing function, keeping only the leading 
terms relevant to the inclination evolution of the jth planet, reads 

^2„2 4 

/t.U- ^^ — V 772^ - 7 (1) / \ r 2 2 / \1 

Rj = — — 2^ ^ , ^, «ifc«jfc^3/2Kfc) bi +Pj - ^{mk + PjPk)\ , 

k=l,^j ® ^ 

where hiLi^^jk) is a Laplace coefficiently rrij is the mass of the jth planet. 



^3/2 V"J 



ajk 



dkl (i-j if cij > ctfc 
aj/ak if aj < a^ 



^ If we declare ajk — oijk — cij/ak and allow the Laplace coefficient to take Ujk > 1 as an argument, 
then separating the case aj < ak from qj > ak is not necessary. We stick here with the textbook convention, 
however. 
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and 

( 1 if a,- > Ofc 

[ ttj/ak it % < Ofc. 

Eqns. (l)-(2) yield two coupled systems of first-order differential equations: (pi, P2, Ps, Pa)'^ = 
A{qi, q2, qs, q^)'^ and (g'l, ^2, g'a, ^4)^ = — ^(pi,P2,P3,P4)"^5 where A is a 4x4 matrix of constant 
coefficients depending on the masses and semi-major axes of the planets. These equations 
describe coupled harmonic oscillators; their solution is 

4 
Qj = ^ Ijk cos(/fct + 7fc) (3) 

k=l 
4 

Pj = ^ Ijk sin(/fct + 7fc) , (4) 

k=l 

where the frequencies fk are the eigenvalues of A. The elements Ijk of the eigenvectors of A, 
and the phases 7^, are fitted to the initial inclinations and ascending nodes of the planets. 
We take initial conditions and planetary data from the NASA JPL Horizons database for JD 
= 2451544.5 (Jan 1 2000; t = in equations (3) and (4)). Table [T] lists the resultant values 
for fk, Ijk, and 7^. We refer to this solution as the BvW solution, even though Brouwer and 
van Woerkom (1950) included all eight planets. 

We now turn to our main concern, the motion of a KBO of negligible mass with semi- 
major axis a > ttj. Its disturbing function is 

2 2 4 

^ = -^ E ^^^3'/2(«./«)b' +P'- 2(9?, + PP,)] , 

where unsubscripted variables refer to the KBO. Because all masses and semi- major axes 
are fixed for this secular problem, Lagrange's equations of motion are of the form 

4 
= -CoP + Yl ^iP^ 

4 

P= — 2^~ = +^0^ ^ Yl ^i^i ' 

i=i 

where the c's are constants. Substituting each equation into the time derivative of the other, 

we find 

4 

g = -c^g + ^ 9 [coqj (t) + Pj (t)] (5) 



1 dR 


no? dp 


1 dR 
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P = -cIp + Y^ Cj \cqVj (t) - (jj (t)] . (6) 

Eqns. ([3])^® describe harmonic oscillators of natural frequency Cq, forced by the planetary 
terms in the sums. The motion is composed of a forced oscillation and a free oscillation: 

q = gforced + gfree = forced + ^free COs(/t + 7) (7) 

P = Pforccd + Pfrec = Pforced + «frec sin(/t + 7) , (8) 

where ifree and 7 are constants determined by initial conditions, and the free precession 
frequency 

/ = Co = 5^ C, = -^ 5^ ^7^3'/2(«j/«) • 

The functions gforced and Pforced depend only on planetary parameters and the KBO semi- 
major axis a: 

4 

^forced = - ^ , _ , COs(/fct + 7fc) (9) 

fc=l ^ ■''' 
4 

Pforced = - J^ ~n~r ^^^^^''^ ^ '^^^ ^^^^ 

k=l •' •'^ 



where 
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The inclination vector i of a KBO is the vector sum of a forced inclination vector 
^forced = (gforced,Pforced) and a free inclination vector ifree = (gfree,Pfree)- Throughout this pa- 
per, we refer to a forced inclination ^forced = (^Lced + Pforced)^''^ ^ forced node ^forced = 
arctan(pforcod/gforcod), a free inchnation z^ee = (gLe + PheeY^"^^ and a free node fifree = 
arctan(pfi.cc/5'frcc)- Thus free and forced inclinations refer to magnitudes (not vectors). The 
free inclination is constant; it is the undamped amplitude of the free oscillation. 

Figure [1] shows the evolution of iforced at various semi-major axes in the Kuiper belt. 
The open diamonds, located at the endpoints of iforced, niark what we call the forced poles. 
At fixed time, the forced poles all lie along a line intersecting the invariable pole, denoted I. 
To the extent that the k = 4 mode, driven mainly by Neptune, dominates, the forced poles 
rotate clockwise (regress) about the invariable pole at a single frequency (f^), maintaining 
constant distance to I. Note that the forced poles at a > 40.5AU lie diametrically opposite 
to those at a < 40.5AU, and that as a approaches 40.5AU from either above or below, the 
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Fig. 1. — Locations of forced poles (open diamonds) at various semi-major axes, at t = (JD 
= 2451544.5) and t = 0.3 Myr. Coordinates are referenced to the ecliptic and mean equinox 
of J2000 (the origin marks the location of the J2000 ecliptic pole). The invariable pole is 
marked by "I." Note how the forced poles at a > 40.5 AU lie diametrically opposite to those 
at a < 40.5 AU, reflecting forcing by the z/ig resonance at a = 40.5 AU. As a approaches 
40.5 AU, the forced pole diverges from the invariable pole. As a increases beyond 40.5 AU, 
the forced pole approaches the invariable pole. We show in this paper that the forced poles 
point normal to the Kuiper belt plane. 
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separation between the forced pole and I increases. These latter two properties reflect the 
z^ig secular resonance at a = 40.5AU, where the denominators / — /4 of equations fl9l)- flT0l) 
vanish: the forced response becomes infinite in magnitude at resonance, and changes phase 
by 180° across resonance (the sign oi f — f\ switches across a = 40.5AU). 

Figure [2] illustrates how free nodes Qf^ee phase mix and how such phase mixing helps 
determine the KBP within the BvW theory. At t = 0, one hundred KBOs having semi-major 
axes within 0.5 AU of 43 AU are set down with orbit normals approximately aligned about an 
arbitrary direction (the orbit poles are actually distributed in a small box in p-q space). Over 
tens of Myrs, the orbit poles of the particles drift away from one another: the small dispersion 
in semi-major axis produces a small dispersion in the free precession frequency / (the rate 
at which the free inclination vector rotates about the forced pole)o While the free nodes 
distribute themselves over all phases, the free inclinations Zfree remain fixed. Thus, at late 
times, the collection of free inclination vectors ifree are distributed axisymmetrically about 
the mean forced pole at a ~ 43 AU. The axisymmetry arises practically independently of how 
the particles' orbit normals are initially distributed; the only requirements are that many 
particles share the same ifree (so that such particles, when phase mixed, trace a full circle) 
and that there exists a small but non-zero dispersion in semi-major axis (so that there exists 
a non-zero spread in free precession frequencies, enabling phase mixing). Note in Figure [2] 
how the invariable pole does not coincide with the center of the circular distribution. 

We conclude that according to the BvW theory, the mean orbit normal of a phase-mixed 
group of particles having approximately the same semi- major axis is given by the forced pole 
corresponding to that semi-major axis. That local forced pole varies with time (Figured]), 
but the particles always encircle it (Figure [2]). The KBP is warped because the forced pole 
changes direction with KBO semi-major axis (Figure [1]). 

Note finally that the ability of test particles to keep their free inclinations constant 
relative to a time-variable forced pole should not be confused with adiabatic invariance. 
Generally the frequencies / and fk are not cleanly separated. The constant if^ee simply 
reflects the undamped amplitude of free oscillation, and is set by initial conditions. In the 
next section, we use numerical simulations to test the constancy of free inclination. 



^In the linear theory, precession frequencies / and fk do not depend on eccentricity or inchnation, but in 
higher order theories they do. 
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Fig. 2. — Phase mixing and axisynimetry of KBOs about the forced pole. In panel (a), the 
orbit poles of 100 test particles (solid circles), having semi-major axes randomly distributed 
between 42.5 and 43.5 AU, are distributed on a grid in p-q space. The forced pole at a = 43 
AU (open diamond) and the invariable pole ("+") are indicated. As time elapses from panels 
(b) through (d), the poles of the test particles precess about the forced pole, which itself 
precesses about the invariable pole. Precession rates differ from one particle to the next 
according to their semi-major axes, while the distance from each pole to the forced pole — 
the free inclination (see panel c) — remains constant. Consequently, after a few tens of Myrs, 
orbit poles tend to phase mix into an annulus centered on the forced pole (not the invariable 
pole). The evolution shown here was computed using the analytic BvW solution, not the 
numerical integration. 
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3. NUMERICAL INTEGRATIONS 

3.1. Initial Conditions 

We calculate the evolution of the KBP at three semi-major axes: a = 38, 43, and 
44 AU. These are chosen to lie away from strong mean-motion resonances (e.g., the 3:2 
resonance resides at 39.5 AU) and outside the a = 40-42 AU region of instability carved by 
the z/18, 1^17 and ug secular resonances (see, e.g., Chiang et al. 2007). At each a we lay down 
Ni X Nq test particles whose initial free inclination vectors are distributed axisymmetrically 
about the local forced pole. That is, each particle's initial p = Pforced(«) + Pfree and initial 
q = qiorced{a) + qtvee, where kree = (pLe + ^lee)^'^^ ^akcs 1 of A^i = 4 values (0.01, 0.03, 
0.1, 0.3 rad) and fifrce = arctan(pfi.ee/(?free) takes 1 of Nq = 20 values distributed uniformly 
between and 2tt. This set-up permits us to directly test the BvW theory, which predicts 
that all Nq particles corresponding to a given a and given initial Zfree should keep the same 
Zfree for all time: a circle of points in p-q space should continue to trace the same-sized 
circle. Coordinates Pforced and ^forced for the initial forced poles are computed using the BvW 
solution of ^ for t = 0. 

Initial osculating eccentricities are zero and initial mean anomalies are chosen randomly 
between and 27r. The four giant planets are included in the integration, with initial 
conditions taken from the JPL Horizons database for JD = 2451544.5 (t = in the BvW 
theory). The integration is performed with the swift_rmvs3 code, written by Levison & 
Duncan (1994) and based on the algorithm developed by Wisdom & Holman (1991). The 
duration of the integration is 4 Gyr and the timestep is 400 days (about 1/11 the orbital 
period of Jupiter). We work in a heliocentric coordinate system, the better to compare with 
the linear secular theory which uses heliocentric elements. 



3.2. Results 

According to the BvW theory, each set of Nq = 20 particles having the same initial 
a and initial i^ee should trace the perimeter of a single circle in p-q space at any given 
time, with the center of the circle yielding the local normal to the KBP. Figure [3] tests 
this prediction; the panels display the p-q positions of particles initially having a = 43 AU 
and ifree = 0.1 rad, sampled at four different times. Most particles at a given time do lie 
approximately on one circle, although there are outliers (see crosses in panels b, c, and d). 
The outliers represent particles whose inclinations and eccentricities grow to large values. 
Many of these particles suffer close encounters with Neptune, whereupon they are removed 
from the integration. We did not identify the cause of the instability, but probably the 
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various high-order mean-motion resonances in the vicinity (Nesvorny & Roig 2001) are to 
blame. Of the 20 particles shown at t = in panel a, only 12 survive to t = 4 Gyr in panel 
d. Those that survive have semi-major axes that remain constant to within ±0.5 AU. 

At each t, we fit a circle (in a least squares sense) to those particles known to survive 
the entire length of the integration. The fit is substantially improved by also discarding 
particles whose instantaneous eccentricities exceed some value Ccut, as we find that particles 
on eccentric orbits tend also to be outliers in p-q space. The value of Ccut is reduced from one 
until the fit parameters cease to change significantly. Note that a particle that is discarded 
from the fit by the Ccut-criterion at one time can be restored to the fit at a later time (if its 
eccentricity falls below Ccut at that later time) . Table [2] lists the values of Ccut chosen for the 
various combinations of initial a and initial Zfree- 

The circles so fitted are overlaid in Figure O A typical fit is excellent and the local 
normal of the KBP (center of the fitted circle; solid diamond) is confidently identified. 
Figures HI [5], and [6] report analogous results for other choices of initial a and initial if^ee- 
All 20 particles having a = 43 AU and ifree = 0.3 rad are stable for 4 Gyr (Figure H]). By 
contrast, at a = 38 AU, where various high-order mean-motion resonances are known to 
cause instability (Nesvorny & Roig 2001), only 2 out of 20 particles having ^free = 0.3 rad 
remain at the end of 4 Gyr (Figure [6]). Fitting a unique circle to 2 points is impossible. And 
as is clear from Figure [6l even if we were to try fitting circles at earlier times when more 
particles are present, such fits would be poor. Only certain combinations of a-ifree permit 
determination of the KBP, as documented in our Table [21 

Panels d of Figures [7] and [8] plot the radii of the fitted circles versus time, for two choices 
of initial a-ifree- The BvW linear theory predicts that these radii should be constant. In fact 
they are nearly so, varying by at most one part out of six. 

How can we further test the BvW theory when we know that it fails to predict the 
orbits of the planets on Gyr timescales? We compute instead a semi-analytic, BvW-based 
solution as follows. At each time in the numerical integration, we output the inclinations 
and ascending nodes of the giant planets and use these to recompute the eigenvectors — and 
thus the forced poles — of the linear theory. Thus we obtain a prediction of where the KBP 
should reside according to the linear theory at a given instant, using the simulation results for 
the planets at that instant to supply the integration constants. This "continuously updated 
BvW solution" for the local normal (forced inclination vector) is shown as an open diamond 
in Figures [SHBl Its location tracks that of the numerically fitted pole (solid diamond) well — 
much better than does the invariable pole (upright cross); see especially Figures [3] and O 
Figures [7| and [8] also demonstrate that the continuously updated BvWP hews closely to the 
numerically fitted KBP. 
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Fig. 3. — Tracking the Kuiper belt plane by numerical integration. In panel (a), we set up 
Nq = 20 test particles having initial semi-major axes of ax = 43 AU and having orbit poles 
distributed in a cone of half- width i^^cc = 0.1 rad centered on the local forced pole given 
by the BvW solution. In panels (b) through (d), the orbit poles evolve according to our 
numerical integration. Solid circles denote test particles that both survive the entire 4 Gyr 
duration of the integration and have osculating eccentricities less than Ccut = 0.08; these 
are fitted to a circle, whose center yields the local normal to the Kuiper belt plane (solid 
diamond). The "X" symbols denote test particles that do not satisfy these requirements. 
The Kuiper belt pole so obtained follows closely that predicted by the continuously updated 
BvW solution (open diamond), and does not point along the invariable pole ("+"). 
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Fig. 4. — Same as Figure [3l except that initial ax = 43 AU, initial Zfree = 0.3 rad, and no 
Ccut criterion is applied. The numerically fitted Kuiper belt pole (solid diamond) tracks the 
continuously updated BvW pole (open diamond) well; both precess about the invariable pole 
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Fig. 5. — Same as Figure [3l except that initial ok = 38 AU, initial ifree = 0.1 rad, and no Ccut 
criterion is applied. Here again the numerically determined Kuiper belt pole (solid diamond) 
is well described by the continuously updated BvW pole (open diamond), not the invariable 
pole ("+")• 
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Fig. 6. — Same as Figure [HI except that initial ok = 38 AU, initial if^ee = 0.3 rad, and no Ccut 
criterion is applied. Only two test particles survive the 4 Gyr duration of the integration, 
rendering determination of the Kuiper belt pole impossible. 
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To summarize the results of our numerical simulations: Linear secular theories like BvW 
correctly predict how the warped KBP evolves with time, provided the parameters of those 
analytic theories are continuously updated using either observations or numerical simulations 
of the giant planets' orbits. The KBP today is accurately predicted by the updated BvW 
solution. The KBP does not, in general, coincide with the invariable plane, except at infinite 
distance from the planets. 



4. THEORY VERSUS OBSERVATION 

We have shown by numerical simulations in §3] that the Classical KBP is given by the 
BvWP, i.e., by linear secular theory. Here we assess whether observations of KBOs bear out 
this result, by locating the actual poles of the KBP near 38 AU and 43 AU. These semi- 
major axes lie to either side of the uis resonance; theory predicts that the corresponding 
poles should lie to either side of the invariable pole, with all three orbit normals lying in 
one plane (see Figured]). Our dataset consists of KBOs listed on the Minor Planet Center 
website on Jan 22 2008 whose (a) astrometric arcs extend longer than 50 days (many objects 
in our sample have much longer arcs), (b) eccentricities are less than 0.1, (c) inclinations are 
less than 10 degrees, and (d) are classified by the Deep Ecliptic Survey (DES) as "Classical" 
in at least two out of their three orbital integrations (see Elliot et al. 2005 for a description 
of the classification scheme). In the vast majority of cases, all three DES integrations yield 
a classification of "Classical." Moreover, the typical 3a uncertainties in semi-major axes are 
smaller than ~0.1 AU. We assemble two samples, one for which 38.09 AU < a < 39.10AU 
(the "38 AU" sample, containing 10 objects) and another for which 42.49 AU < a < 43.50 AU 
(the "43 AU" sample, containing 80 objects). Object designations are given in Table O 

Figure [9] plots the observed {q,p) positions for the two samples, and compares with 
theory. Theory predicts that the average {q,p) measured for each sample should equal 
(ftorced^Pforced) Calculated for the average semi-major axis of the sample. The good news is 
that the data are consistent with this prediction; the differences between the observed poles 
and the forced poles are less than 2a for both samples (the error bars in Figure [9] are ±3cr, 
where a is the standard deviation of the mean). The bad news is that the observed poles are 
also consistent with the invariable pole, at a similar confidence level. While it is encouraging 
for the theory that the observed pole at 38 AU indeed lies to the left (toward smaller q) 
of the invariable pole, and that the observed pole at 43 AU lies to the right, there are not 
enough observations to make more precise statements and to rule out the hypothesis that 
the KBP equals the IP with greater than 3a confidence. 

Note that by selecting our sample to have orbital inclinations less than 10 degrees with 
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Fig. 7. — Results for test particles having initial ok = 43 AU and initial ifree = 0.1 rad, 
demonstrating that the numerically obtained KBP follows the continuously updated BvWP, 
not the IP. (a) Inclination of the KBP (solid circles), BvWP (open circles), and IP (line), all 
relative to the ecliptic, (b) Mutual inclination between the KBP and BvWP (solid diamonds), 
the BvWP and the IP (open diamonds), and the KBP and the IP (solid line), (c) Longitude 
of ascending node of the KBP (solid circles), BvWP (open circles), and IP (line), all relative 
to the ecliptic and mean equinox of J2000. (d) Numerically obtained free inclination (solid 
circles; these are the radii of the fitted circles in Figure [3]), compared against the initial ifree 
(line). The BvW theory predicts that the free inclination should be constant; while it is 
nearly so in our numerical integration, its mean value is offset by 0.6 deg (10%) relative to 
its initial value. Data are sampled every 10^ yr and do not resolve the precession of the KBP 
occurring with a ~10^-yr period (see Figure [TOj) . 
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Fig. 8. — Same as Figure [71 except for initial qk = 38 AU and initial if^ee = 0.1 rad. 



-18- 



CD 

CD 
"O 



10 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

• 


5 


• 

• 

- o - 

• o 

t • 




•; • . . 




• O: O . 

cP o'i • •• 




A* ^..TV: • • • 




^.''U^" . 




•0.*./ . 





■■'■■■ °* *. (3 • • 




• .. • 




o 




• 






. += invarloble pole 


• 


-5 


- ^~ 'ovg,38^~ ovg,4j 


• 




- ♦= 'forced (43.0AU) 


- 




- 0= 'forced (38.6AU) 


• 




. • (o=43.0±0.5AU) 


• 




. o (o=38.6±0.5AU) 


_ 


10 


. . . . 1 . 


. , 1 .... 1 ... , 



-10 



-5 

q (degree) 



10 



Fig. 9. — Pole positions of observed Classical KBOs having semi-major axes near 38 AU 
(open circles) and 43 AU (filled circles), in g = icosQ,p = isinQ space, referenced to the 
J2000 ecliptic. According to theory (||2]-fj3]), the average pole position at a given semi-major 
axis (triangles with error bars) should match the forced pole positions (diamonds) calculated 
from BvW. They do to within 2a (error bars are ±3cr, where a is the standard deviation of 
the mean). Unfortunately, the alternate hypothesis that the average pole positions are given 
by the invariable pole (bold cross) cannot be ruled out with greater confidence. 
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respect to the ecliptic plane, we bias our measurement of the average pole position towards 
the ecliptic pole. This systematic error is probably still smaller, however, than our random 
error. For example, the observed pole at 38 AU actually lies further from the ecliptic pole 
than does the theoretically expected forced pole. See Elliot et al. (2005) for ways of reducing 
this systematic bias. 

We also performed a two-dimensional Kolmogorov-Smirnov test to see whether the (g, p) 
distributions for the two samples differ (they should). The probability that they do not is 
4.9% — small enough to be suggestive of a real difference, but in our judgement too large to 
be conclusive. 



5. SUMMARY AND DISCUSSION 

Classical Kuiper belt objects trace a sheet that warps and precesses in response to the 
planets. The current location and shape of this sheet — the "plane" of the Classical Kuiper 
belt — can be computed using linear secular theory, with the observed masses and current 
orbits of the planets as input parameters (Brown & Pan 2004). The local normal to the 
plane is given by the theory's forced inclination vector, which is specific to every semi-major 
axis. At infinite distance from the planets, the plane coincides with the invariable plane. The 
deviations of the Kuiper belt plane away from the invariable plane, while generally non-zero, 
are typically small: less than 3 degrees outside the secularly unstable gap at a ~ 40.5 ±1 AU 
(inside the gap no KBOs have been observed, as expected). As the semi- major axis varies 
from < 40.5 AU to > 40.5 AU, the ascending node of the Kuiper belt plane on the invariable 
plane rotates by very nearly 180 degrees, a result of the sign change in the forced response 
across the i^ig resonance. 

These conclusions are supported by our numerical integrations of giant planet and test 
particle orbits lasting 4 Gyr. These integrations show that a Kuiper belt object maintains 
a nearly fixed orbital inclination with respect to the time-variable Kuiper belt plane. This 
accords with linear theory, in which the free inclination represents the undamped, constant 
amplitude of a test particle's free oscillation. It may seem surprising that the linear theory 
is vindicated in this regard while it cannot accurately predict planetary motions and hence 
the future location of the Kuiper belt plane. But the inaccuracies accrue only slowly — 
Figure [10] shows that it takes several precession periods for the linear theory to diverge from 
numerical simulation in predicting the location of the plane, and the analytic and numerical 
solutions are always qualitatively similar. Referring back to the equations of motion ([5])-([6]), 
we see that if we take the planetary forcing terms Qj (t) and pj (t) to be given by the more 
realistic numerical integration — instead of the inaccurate but qualitatively similar analytic 
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solution Q-dl]) — then the free component of the motion (the homogeneous solution of the 
differential equation) would still be given by dZ])-®, irrespective of the forced component 
(the inhomogenous solution). Thus an object can keep its free inclination with respect to 
the forced plane fixed, even though the location of the forced plane itself cannot be forecast 
analytically. 

Currently the data on actual KBOs are consistent with, but do not conclusively verify, 
our theoretical finding that the Kuiper belt plane warps by a few degrees to either side of 
the secularly unstable gap. Our analysis of the observations, like that of Elliot et al. (2005), 
suffers from large random errors. The study by Brown & Pan (2004) does not, but at the 
expense of including objects of all dynamical classes and not following variations in the pole 
position with semi-major axis. Quadrupling the sample size of low-i, low-e objects at 38 
AU, where there are currently ten usable objects, can increase our confidence in the reality 
of the warp to greater than 3a. 

Detecting "Planet X" via its influence on the Kuiper belt plane will be substantially 
more challenging. We calculate that a 100-M^ planet with a semi-major axis of 300 AU and 
an orbital inclination with respect to the ecliptic of 10 degrees would shift the forced poles 
in the 43-50 AU region by 0.1 degree. 

We thank Ruth Murray-Clay and Mike Brown for discussions, and Chris Culter for sug- 
gesting this undergraduate research collaboration. An anonymous referee provided numerous 
helpful comments. This work was supported by NSF grant AST-0507805. 
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Table 1. BvW solution" for Jan 1, 2000 (t = 0) 



k = 1 



1 


1.5792 


0.36081 


0.055195 


-0.066653 


2 


1.5792 


-0.89904 


0.044853 


-0.064203 


3 


1.5792 


0.040670 


-1.01200 


0.062630 


4 


1.5792 


0.0045122 


0.11941 


0.67290 



fk 

Ik 



1.87J 



-123.13 
-0.91730 



-14.102 
2.3699 



-3.2984 
-2.7418 



'^The components Ijk of the eigenvectors (in de- 
grees), eigenfrequencies fk (in radians per Myr), and 
phases 7^ (in radians), calculated using data for the 
giant planets on JD = 2451544.5 from NASA JPL 
Horizons. All elements are heliocentric and referred 
to the ecliptic and mean equinox of J2000. 
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Table 2. Numerical Results for Stability in Initial a-ifree Space 



a 


^frcc 


7^ of Survivors 


Cmax 


Ccut 


Permits KBP 


(AU) 


(rad) 


(out of 20) 


of Survivors 




Determination? 


38 


0.01 


17 


0.05 


None Applied 


Yes 




0.03 


16 


0.05 


None Applied 


Yes 




0.10 


8 


0.07 


None Applied 


Yes 




0.30 


2 


0.02 


None Applied 


No 


43 


0.01 


19 


0.12 


0.04 


Yes 




0.03 


20 


0.04 


None Applied 


Yes 




0.10 


12 


0.12 


0.08 


Yes 




0.30 


20 


0.02 


None Applied 


Yes 


44 


0.01 


20 


0.11 


0.03 


Yes 




0.03 


20 


0.11 


0.04 


Yes 




0.10 


20 


0.15 


0.06 


Yes 




0.30 


20 


0.29 


0.08 


Yes 
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Table 3. Samples of Observed KBOs for Locating the KBP 



Sample Objects 



'38 AU" 1998 WV24, 1999 OJ4, 2000 YB2, 82157, 2003 FD128, 2003 QA92, 2003 YL179, 

2003 QQ91, 144897, 119951 
^'43 AU" 19255, 1994 EV3, 1996 TKge, 33001, 1998 WY24, 1998 WX24, 1999 CN153, 

1999 RT214, 1999 ON4, 1999 XY143, 1999 RW214, 1999 CH154, 1999 RU215, 

1999 HVii, 1999 DA, 1999 HJ12, 1999 CW131, 2000 PU29, 2000 PX29, 
134860, 2000 CL105, 2000 ONqt, 2000 PC30, 2000 FS53, 2000 WV12, 

2000 WLi83, 2000 OUgg, 2000 YUi, 88268, 2001 QB298, 2001 QD298, 2001 XR254, 

2001 QO297, 2001 HZ58, 2001 RW143, 2001 OK108, 2001 DBioe, 88267, 

2001 XU254, 2001 FKi85, 2001 OZ108, 2002 CD251, 2002 PX170, 2002 PV170, 

2002 FW36, 2002 WL21, 160256, 2002 VB131, 2002 PY170, 2002 CS154, 2002 PD155, 

2003 SN317, 2003 UT291, 2003 FK127, 2003 QG91, 2003 FA130, 2003 HY56, 

2003 QY90, 2003 TK58, 2003 QF91, 2003 QE91, 2003 QZm, 2003 QL91, 2003 QE112, 
2003 YR179, 2003 QYiii, 2003 QD91, 2003 TL58, 2003 QU90, 2003 YT179, 
2003 YX179, 2003 YS179, 2003 YJ179, 2004 UDio, 2004 DM71, 2005 JZ174, 
2005 GDi87, 2005 JP179, 2005 XUioo, 2006 HA123 



